*     potential at CZ, RTIME for areasink with corners CZ1 CZ2 CZ3 CZ4
*     in counterclockwise order with timely and spacially constant
*     extraction for RTIME>0 and 0 for RTIME<0) with transient constant
*     RP (=porosity/(hydaulic conductivity*average head)
      DOUBLE PRECISION
     &   FUNCTION RFSAP1(CZ, RTIME, CZ1, CZ2, CZ3, CZ4, RP)
      IMPLICIT NONE

      COMPLEX*16 CZ, CZ1, CZ2, CZ3, CZ4
      DOUBLE PRECISION RTIME, RP

      DOUBLE PRECISION RPI, RSQRPI, RTOLRL
      PARAMETER (RPI=3.1415926535)
      PARAMETER (RSQRPI=1.77245385090551)
      PARAMETER (RTOLRL=.0001)

      COMPLEX*16       CZL(4), CZN
      DOUBLE PRECISION RHL(4), DFUN, RPO4T, RSQPO4, RTOL,
     &                 RLIM, RX, RY, RX1, RX2, RX1XYY, RX2XYY,
     &                 RFIEER, RFERFC, RFE1, D1, D2, D3, D4,
     &                 DMAX1, CDABS
      INTEGER          I
      LOGICAL          LIN, LFSAIN

      DFUN=0.
      RPO4T =.25*RP/RTIME
      RSQPO4=SQRT(RPO4T)
*     half length of the sides
      RHL(1)=.5*CDABS(CZ2-CZ1)
      RHL(2)=.5*CDABS(CZ3-CZ2)
      RHL(3)=.5*CDABS(CZ4-CZ3)
      RHL(4)=.5*CDABS(CZ1-CZ4)
*     tolerance relative to the longest side
      RTOL=RTOLRL*2.*DMAX1(RHL(1), RHL(2), RHL(3), RHL(4))
*     use CZL(i) for temporary storage of the corner i
      CZL(1)=CZ1
      CZL(2)=CZ2
      CZL(3)=CZ3
      CZL(4)=CZ4
*     get information about CZ being inside or outside and obtain
*     new coordinates CZN (RTOL) away from the sides if CZ is too
*     close to a side
      LIN=LFSAIN(CZ, CZL, RHL, RTOL, CZN)
*     use CZL(i) for local coordinates with respect to side i
      CZL(1)=(CZN-.5*(CZ2+CZ1))*2.*RHL(1) / (CZ2-CZ1)
      CZL(2)=(CZN-.5*(CZ3+CZ2))*2.*RHL(2) / (CZ3-CZ2)
      CZL(3)=(CZN-.5*(CZ4+CZ3))*2.*RHL(3) / (CZ4-CZ3)
      CZL(4)=(CZN-.5*(CZ1+CZ4))*2.*RHL(4) / (CZ1-CZ4)
      DO I=1,4
         RX1=-RHL(I)
         RX2=RHL(I)
         RX=REAL(CZL(I))
         RY=AIMAG(CZL(I))
         IF (.NOT.(ABS(RX).GT.RX2.AND.ABS(RY).LT.RTOL)) THEN
            RLIM=RSQPO4*RY
            RX2XYY=(RX-RX2)*(RX-RX2)+RY*RY
            RX1XYY=(RX-RX1)*(RX-RX1)+RY*RY
            D1= .5/RSQRPI*(.5*RP*RY*RY+RTIME)*
     &         RFIEER( RLIM, (RX-RX2)/RY, (RX-RX1)/RY )
            D2= -.25/RSQRPI*RY*SQRT(RTIME*RP)*EXP(-RPO4T*RY*RY)*
     &         ( RFERFC(RSQPO4*(RX-RX2))-RFERFC(RSQPO4*(RX-RX1)) )
            D3=  .125*(RX-RX2)*RY*RP/RPI * RFE1(RPO4T*RX2XYY)
            D4= -.125*(RX-RX1)*RY*RP/RPI * RFE1(RPO4T*RX1XYY)
            DFUN=DFUN+D1+D2+D3+D4
         END IF
      END DO
      IF (LIN) DFUN=DFUN-RTIME
      RFSAP1=DFUN/RP
      RETURN
      END
